home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / GAUSSJ.DEM < prev    next >
Text File  |  1991-04-29  |  2KB  |  86 lines
  1. PROGRAM d2r1 (input,output,dfile);
  2. (* driver program for subroutine GAUSSJ *)
  3. LABEL 10,99;
  4. CONST
  5.    np=20;
  6.    mp=20;
  7. TYPE
  8.    glnpbynp = ARRAY [1..np,1..np] OF real;
  9.    glnpbymp = ARRAY [1..np,1..mp] OF real;
  10.    glnp = ARRAY [1..np] of integer;
  11. VAR
  12.    j,k,l,m,n : integer;
  13.    a,ai,u : glnpbynp;
  14.    b,x,t : glnpbymp;
  15.    dfile : text;
  16.  
  17. (*$I MODFILE.PAS *)
  18. (*$I GAUSSJ.PAS *)
  19.  
  20. BEGIN
  21.    glopen(dfile,'matrx1.dat');
  22. 10:   readln(dfile);
  23.    readln(dfile);
  24.    readln(dfile,n,m);
  25.    readln(dfile);
  26.    FOR k := 1 to n DO BEGIN
  27.       FOR l := 1 to n-1 DO read(dfile,a[k,l]);
  28.       readln(dfile,a[k,n])
  29.    END;
  30.    readln(dfile);
  31.    FOR l := 1 to m DO BEGIN
  32.       FOR k := 1 to n-1 DO read(dfile,b[k,l]);
  33.       readln(dfile,b[n,l])
  34.    END;
  35. (* save matrices for later testing of results *)
  36.    FOR l := 1 to n DO BEGIN
  37.       FOR k := 1 to n DO BEGIN
  38.          ai[k,l] := a[k,l]
  39.       END;
  40.       FOR k := 1 to m DO BEGIN
  41.          x[l,k] := b[l,k]
  42.       END
  43.    END;
  44. (* invert matrix *)
  45.    gaussj(ai,n,np,x,m,mp);
  46.    writeln;
  47.    writeln('Inverse of matrix a : ');
  48.    FOR k := 1 to n DO BEGIN
  49.       FOR l := 1 to n-1 DO write(ai[k,l]:12:6);
  50.       writeln(ai[k,n]:12:6)
  51.    END;
  52. (* test results *)
  53. (* check inverse *)
  54.    writeln('a times a-inverse (compare with unit matrix)');
  55.    FOR k := 1 to n DO BEGIN
  56.       FOR l := 1 to n DO BEGIN
  57.          u[k,l] := 0.0;
  58.          FOR j := 1 to n DO BEGIN
  59.             u[k,l] := u[k,l]+a[k,j]*ai[j,l]
  60.          END
  61.       END;
  62.       FOR l := 1 to n-1 DO write(u[k,l]:12:6);
  63.       writeln(u[k,n]:12:6)
  64.    END;
  65. (* check vector solutions *)
  66.    writeln;
  67.    writeln('Check the following vectors for equality:');
  68.    writeln('original':20,'matrix*sol''n':15);
  69.    FOR l := 1 to m DO BEGIN
  70.    writeln('vector ',l:2,':');
  71.       FOR k := 1 to n DO BEGIN
  72.          t[k,l] := 0.0;
  73.          FOR j := 1 to n DO BEGIN
  74.             t[k,l] := t[k,l]+a[k,j]*x[j,l]
  75.          END;
  76.       writeln(' ':8,b[k,l]:12:6,t[k,l]:12:6);
  77.       END
  78.    END;
  79.    writeln('***********************************');
  80.    IF eof(dfile) THEN GOTO 99;
  81.    writeln('press RETURN for next problem:');
  82.    readln;
  83.    GOTO 10;
  84. 99:   close(dfile)
  85. END.
  86.